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Abstract 

We address the inverse problem of designing isotropic pairwise particle interaction potentials that lead 
to the formation of a desired lattice when a system of particles is cooled. The design problem is motivated 
by the desire to produce materials with pre-specified structure and properties. We present a heuristic 
computation-free geometric method, as well as a fast and robust trend optimization method that lead to 
the formation of high quality honeycomb lattices. The trend optimization method is particularly successful 
since it is well-suited to efficient optimization of the noisy and expensive objective functions encountered in 
the self-assembly design problem. We also present anisotropic potentials that robustly lead to the formation 
of the kagome lattice, a lattice that has not previously been obtained with isotropic potentials. 
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I. INTRODUCTION 



During the process of self-assembly, randomly distributed components acting under the influence 
of short-range mutual interactions, arrange themselves into a highly ordered final configuration. 
An important feature of self-assembly is that the components arrange themselves into this more 
ordered state without the infiuence of external factors. The ordered arrangement observed at the 
level of the resulting superstructure is predicated on the design of the individual components and 
their local interactions. Understanding the self-assembly process and how local properties of the 
components may be manipulated to infiuence the resulting global ordering, is an active area of 
research spanning a broad range of disciplines including materials science, chemical engineering, 
bioengineering, and nanotechnology. 

In a seminal article, Whitesides and Grzybowski [l| provided a broad definition of self-assembly, 
and addressed promising applications in a wide range of disciplines. They observed that self- 
assembly of cells to form tissue, organs, and ultimately organisms, is fundamental to life, and 
motivates a determined study of the self-assembly process. Subsequent studies and applications of 
self-assembly are as fascinating as they are plentiful. Murr et al 3] identified self-assembly as the 
mechanism by which silicatein monomers combine to form protein fibers in certain marine sponges. 
Zheng et al used shape recognition and selective binding to induce self-assembly and packaging 
of integrated semiconductor microsystem devices by agitating the components in an aqueous en- 



vironment Q], and Stauth et al [J] have manufactured field-effect transistors via self-assembly of 
micrometer scale components using similar techniques. In ^, self-assembly is shown to enhance 
synergistic group transport of autonomous robots. Jakab et al recognized that although biolog- 
ical systems are genetically controlled, the formation of biological superstructures is ultimately 
governed by physical interactions, and demonstrated the formation of prescribed shapes using 
self-assembling multicellular systems 0|. 

Typically, studies of self-assembly examine the ordered superstructures that arise from a given 
fixed interaction potential. For example, Manoharan et al used optical and electron microscopy 
to identify the range of structures produced by self-assembly of colloidal microspheres as fluid is 
removed from the emulsion droplets containing the spheres Q]- Maksymovych et al investigated 
the formation and reactivity of linear chains of dimethyldisulflde molecules on a gold surface [s]. 
Engel (01, see also [ll|, [l^) have observed the formation of complex crystals and 

quasicrystals arising from a simple double-well interaction potential. 

Interest in the fabrication of nanomaterials and photonic crystals with desired material proper- 
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21[ motivates the inverse problem: design the shortrange interactions 



in order to induce self-assembly of the components into a desired lattice structure. Laboratory tech- 



2M, and 



niques now available allow for modification and tuning of particle interaction potentials 
hence experimentalists are achieving ever-increasing control over the local interactions that influ- 
ence the global properties of the material formed via self-assembly. These methods typically use 
colloidal suspensions and optical forcing to alter the chemical environment and screening properties 
of the solution in which the assembly occurs. 

In this paper, we specifically consider the problem of designing short-range pairwise interaction 
potentials between particles on a planar surface to induce self-assembly of a desired lattice. The 
main results to be presented are new methods — a heuristic geometric method as well as a robust 
trend optimization method for the design of isotropic interaction potentials that lead to high quality 
honeycomb lattices as the system of particles is cooled. The geometric method is also extended to 
the case of anisotropic potentials which allows for the formation of more exotic kagome lattices. 
Another contribution of this work is the development of tools for objective assessment of quality 
of lattices, that mimic intuitive human perception of lattice quality. 

Rechtsman et al in 



23| have already demonstrated computational methods for finding solutions 
to the inverse self-assembly problem. They used a simulated annealing optimization procedure to 
find potentials that lead_to the self-assembly of particles into square and honeycomb lattices. To 
be sure, the intent of 



23( 1 was to demonstrate that the inverse problem of potential design for the 
purpose of inducing the formation of a target lattice can be solved in practice, and to carefully 
verify that the potentials they proposed do indeed lead to the target lattices through Monte Carlo 
simulation. Hence, the computational effort required to obtain the potentials was only a marginal 
consideration in their work. Presently, we consider the straightforward simulated annealing method 
of [23] as a baseline method with which we may compare the new methods presented here. When 



compared with the baseline simulated annealing method, the optimization procedure described in 
this paper leads to a hundredfold speed-up in the generation of potentials, as well as the formation 
of higher quality lattices. Furthermore, the procedure for finding potentials is more robust, and the 
resultant potentials form the target lattices more robustly with respect to variations in the initial 
conditions of the particles. As will be demonstrated, the chief reason for the marked speed-up 
over the simulated annealing method is the facility of the trend optimization method to optimize 
objective functions that are both noisy and expensive to evaluate. 

The organization of the paper is as follows. In Section [Til we precisely define the self-assembly 
problem, and summarize the method for generating potentials devised previously by [3] that will 
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serve as a baseline method for purposes of comparison. In Section Hill we establish objective metrics 
for measuring lattice quality so that reasonable comparisons between the methods can be made. 

We approach the self-assembly problem by framing it as an optimization problem in which 
the desired potential optimizes a suitably chosen objective function. This approach requires that 
we choose both an optimization scheme and an objective function to be optimized. Section IIVI 
describes three objective functions that will be used, while Section |V] describes in detail the trend 
optimization method. A discussion on the hierarchical nature of the trend method is also provided. 
We proceed in Section IVII to list the five solution methods to be compared and describe their 
implementation. The final comparison of the methods is presented graphically in the plots of 
Section IVIII In Section IVIIIl we present an extension of the geometric method to anisotropic 
potentials that lead to self-assembly of the kagome lattice. 

All molecular dynamics simulations performed during the design and testing of the potentials 
were executed on the CITerra high performance computing cluster housed in the Division of Geo- 
physical and Planetary Sciences at Caltech using the LAMMPS software package Q from Sandia 
Laboratories. 
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Kevrekidis, Alison Marsden, Matthew West, Jose Miguel Pasini, and Igor Mezic for their interest 
and helpful comments. This work was supported in part by DARPA DSO under AFOSR contract 
FA9550-07-C-0024. 



II. THE SELF ASSEMBLY PROBLEM 



The availability of laboratory methods to tune interaction potentials between components, and 
hence influence the structure of the resulting self-assembled configurations, motivates the use of 
self-assembly to produce materials with desired structural properties. The specific self-assembly 
problem addressed in this paper is defined as follows: 

Definition: The Self- Assembly Problem 

Design a radially symmetric pairwise interaction potential, Vhc('^)i so that when a 
system of particles interacting with each other in the plane through this potential is 
cooled, the particles form a honeycomb lattice. 



The purpose of the present paper is to compare methods for generating potentials that solve the 
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self-assembly problem. The methods for generating potentials are compared using three criteria: 



CI. The computational effort required by the method to produce the potential; 
C2. The quality of the lattices formed by the potential; 

C3. The robustness of the quality of the lattices formed to variations in the initial conditions of 
the particles. 

Certainly, for a given method we expect to see trade-offs between these criteria. For example, 
a faster method may lead to a potential that produces lattices of poorer quality. 



Laboratory techniques for tuning interaction potentials motivated Rechtsman et al [23f| to con- 
sider physically realizable potentials as basis functions for the desired potential, where the basis 
functions contain parameters that allow for tuning the shape of the total potential obtained from 
their sum. For the case when the desired final configuration of particles is the honeycomb lattice, 



23] proposed an interaction potential consisting of the sum of a Lennard-Jones potential, an ex- 
ponentially decaying potential, and a Gaussian-shaped potential, parameterized in the following 
way: 

FHc(r; ao, ai, a2, as) = ^ - ^ + ai e""^^^ - ^Ae-^^^^-'^'^^", (1) 

where ao, ai, a2, and as are four free parameters that can be tuned to adjust the shape of the 
potential. 

A solution (there are many) to the self-assembly problem has been provided by [23(]: namely, 



choosing ao = 5.89, ai = 17.9, a2 = 2.49, and as = 1.823 in the expression for Vhc- A sample lattice 
obtained when cooling a system of particles using these parameters for the interaction potential 



is shown in Figure 1(c) The honeycomb lattice is the dominant structure in the lattice, although 
there are still visible defects that arise due to the finite duration of the cooling schedule. The 
cooling simulation used to obtain this lattice, as well as all other simulations referred to in this 
paper, were performed using periodic boundary conditions. 

A reasonable question to ask at this point is: "How difficult is the self-assembly problem?" 
Experience shows that although easy to state, the self-assembly problem is difficult to solve in that 
solutions that lead to the honeycomb lattice are difficult to find and possibly very fragile. For 



instance, adding small perturbations to the parameters in 



231 ] for the honeycomb potential leads 



to the configuration in Figure 2(b) in which the honeycomb structure is less pronounced. The 
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(a) A perfect honeycomb lattice, (b) The potential proposed in [23]. (c) A sample lattice obtained 

using the potential shown in (b). 

FIG. 1: The self-assembly problem entails finding an interaction potential that leads to the formation of a honeycomb 
lattice (a). [2^ used a simulated annealing optimization procedure to generate a potential (b). A sample lattice 
formed using this potential is shown in (c) and exhibits defects that result from the finite duration of the cooling 
simulation. 




(a)ao = 5.89, ai = 17.9, 
a2 = 2.49, as = 1.823. 





(b) ao = 5.89, ai = 17.9, 
a2 = 2.49, as = 1.89. 




(c) ao = 5.0, ai = 17.0, 
a2 — 2.0, as = 1.5. 



(d) ao = 6.0, ai = 18.0, 
a2 = 3.0, aa = 2.0. 



FIG. 2: The lattices shown here are obtained by slowly cooling 484 particles using the interaction potential 



Vhc('"; fflo, ai, 02, as) provided in Equation [T] with the parameter values indicated. Figure 2(a) uses the parameters 



proposed by [23l |. while Figure 2(b) indicates how fragile the honeycomb lattice is with respect to slight changes in 



the parameters. Figures 2(c) and 2(d) indicate the range of structures obtained for generic values of the parameters. 



amorphous configurations shown in Figures 2(c) and |2(d)] indicate the range of structures that are 
possible for generic values of the parameters. 

We have identified the honeycomb lattice as the target lattice in the self-assembly problem 
because it the simplest lattice structure that represents a non-trivial case for self-assembly. The 
triangular lattice is easily and robustly assembled by a straightforward Lennard-Jones interaction 
potential. This fact can be readily demonstrated through direct simulation; however, it is interest- 
ing that a rigorous mathematical proof that the triangular lattice is the global energy minimizer 
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and ground state of a particle system with a pairwise Lennard-Jones interaction potential was 
provided only recently by Theil in 



25| . Generating a potential that robustly induces the formation 
of a square lattice in particle simulations is also a simple matter using the geometric method to be 
presented later. Producing the honeycomb lattice using only a radially symmetric potential is more 
troublesome chiefly because regions of the domain tend to form the triangular lattice - a lattice 
that competes strongly with the honeycomb lattice. This is due to the fact that the triangular 
and honeycomb lattices have identical distances to nearest neighbors, and the triangular lattice is 
formed by simply adding a particle to the center of each hexagonal cell in the honeycomb lattice. 

For concreteness, we constrain the parameter space over which we search for parameter values 
in V\-\c{r;ao,ai,a2,a3) such that 

4.0 < ao < 8.0, < ai < 30.0, < aa < 3.8, and 1.25 < as < 2.25 . (2) 

These ranges yield a large class of function shapes over which the search is performed. 

gori thms to find potentials that lead 
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Rechtsman and co-workers developed two computationa _ 
to the self-assembly of particles into a given target lattice 

The first optimization scheme chooses the shape of the potential so that the energy difference 
between the target lattice configuration and the configuration of the competitor lattices is max- 
imized. The optimization is performed while ensuring mechanical stability by allowing only real 
phonon frequencies. This method is purely static in that it seeks to ensure that the target lat- 
tice configuration is energetically the most preferred final state; the method does not incorporate 
information about the dynamics of the particles as they tend towards this final configuration. 

The second optimization scheme considered by Rechtsman et al concentrates on choosing poten- 
tials so as to maximize the stability of the target lattice near its melting point, while also requiring 
stability of the the target lattice with respect to changes in density, and mechanical stability by 
ensuring that the phonon frequencies are real. After placing particles into the target lattice con- 
figuration, a short molecular dynamics simulation is performed at a fixed temperature just below 
the melting temperature of the lattice. The deviation of the the final configuration from the initial 
target lattice is computed using the Lindemann parameter, 



LP := 



^E('— rT- iE('.-n) > (3) 



where N is the number of particles, and rf'^ and are the initial and final positions of particle i, 
respectively. A simulated annealing procedure is used to choose parameters in the expression for 
HhcI^^; flO) oi) 02, 03) that minimize the Lindemann parameter. 



After generating a potential using these methods, Rechtsman et al carefully checked using Monte 
Carlo simulations that particles starting in a random initial configuration do indeed self-assemble 
into the target lattice. The determination as to whether or not the target lattice was formed, 
was made by visual inspection of the final configuration and deciding if the configuration has few 
enough defects to be considered a lattice, as well as checking the long range order by visually 
inspecting simulated Bragg diffraction patterns. 



III. LATTICE QUALITY METRICS 

As described in Section [III one of the criteria used to compare methods for generating interaction 
potentials is the quality of the resulting lattices. In the work of {23], a simple visual check of the 
resulting lattice was used to determine if the desired lattice (with a few defects perhaps) was 
obtained, and since their purpose was to show that the desired lattices can be obtained, this visual 
check was sufficient. The purpose of the present paper is to compare several methods according to 
the lattice quality criterium, thus we must first introduce objective methods for assessing lattice 
quality. 

A prevalent metric for lattice quality is the structure factor, a quantity that assesses long range 
spectral order in the lattice using diffraction patterns. We have found that this quality metric is 
inadequate for our purposes, firstly because the structure factor does not provide a scalar value for 
quality, and secondly because we have observed that long range order is not necessarily strongly 
correlated with visual perception of lattice quality - a lattice with glaring defects may still exhibit 
a high degree of long range order, for example, while a lattice that has weak long range order 
due to a grain boundary may have excellent local lattice structure. Hence, we have developed two 
lattice quality metrics that mimic nearly as possible the visual assessment of lattices. When the 
human eye looks at a lattice and makes a judgement with respect to lattice quality, the emphasis 
is on order within sub-regions of the entire lattice. A lattice that consists of two perfectly formed 
sub-lattices that have a domain wall where they meet will be judged by the eye to be quite well- 
formed. Thus, although long range order is important, a determination of local ordering is crucial 
for assessing lattice quality in a manner similar to the eye. The metrics that we use to quantify 
lattice quality have this local feature. 

The two lattice quality metrics we present here are called the Template Measure and the Defect 
Measure. In both lower value of the metric corresponds to a higher quality lattice. 
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A. Template Measure 



The Template Measure uses a small segment of the target lattice as a template with which 
to locally compare nearby lattice particle positions for each particle in the given lattice. In the 
honeycomb lattice, a suitable template may be one hexagonal cell composed of six particles, or 
points. For each particle in a given lattice configuration, the first point in the template is pinned 
to the particle, and the template is then rotated to find the best fit to other nearby particles. 
As the template is rotated, each point in the template is paired with the nearest particle in the 
given lattice. The angle of rotation of the template that produces the least deviation between the 
template points and lattice particles is considered the best fit. Once this best fit position of the 
template has been found for each particle in the lattice, the Template Measure (TM) is obtained 
by summing the deviation in the positions of the template points and lattice particles from the 
best fit for each particle in the lattice: 

N 



100 . 
TM := > mm 

N ^ Pi 



N 

p=l 



, 2 

''P, extra 



(4) 



i=2 

where the index p ranges over all N particles in the given lattice, 9 is the angle of rotation of the 
template, j.^'*''™p'^*° position of the z*^ point in the template when the template is attached 

to particle p and rotated by angle 9, rf^ is the position of the particle in the given lattice that 
is closest to j.^'^'^^pi^*'^^ g^Yid c is the number of points in the template (c = 6 for the honeycomb 
cell template). Notice that since the first point in the template is pinned to the given lattice 
particle their positions are equal, that is r^^ = j.^'**^™?!'^*^^ g^j^^j hence the sums over the template 
points need not consider this first template point. The extra term, rip^extra; is a count of any extra 
particles in the given lattice that fall inside the hexagonal template, but are not paired with any of 
the template points. In this way, the best fit of the template seeks not only to match the particle 
locations, but also the void within the hexagon. This 'opacity' of the template ensures that defects 
that arise due to the formation of the triangular lattice that has a particle located at the center of 
the hexagon will be penalized. The prepended scaling factor of is not strictly necessary in 

the Template Measure, but is included for convenience so that the resulting lattice quality values 
can be interpreted as a measure of the defectiveness per particle, and have a magnitude in a range 
from zero to roughly 100. 

An illustration of how the Template Measure is implemented in practice is provided in Figure [3l 



In Figure 3(a) , a single hexagonal honeycomb cell template is attached to a particle in a honeycomb 
lattice and rotated until a best fit with the surrounding particles is achieved. Repeatedly attaching 
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the template to each particle in the lattice and realigning as in Figure |3(b)| quickly reveals the 
locations of the defects. 




(a) (b) 

FIG. 3: Illustration of the Template Measure, (a) A template consisting of a single hexagonal cell is pinned to a 
particle in the given lattice, and rotated to the position which minimizes the distance between points in the template 
and the nearest particles in the lattice, (b) Fitting the template cell to each particle in the lattice and rotating to 
find the best fit, quickly reveals the locations of defects. 



B. Defect Measure 

The second lattice quality metric we present is the Defect Measure. The Defect Measure provides 
a weighted count of all the defects in the local neighborhood of each particle in the given lattice. 
The types of defects considered are shown in Figure HI and include displaced, missing, and extra 
particles. Note that all of the possible defects, including global defects, are taken into account by 
the types of defects shown. For example, extended grain boundaries are taken into account by 
contributions to the Defect Measure from locally displaced, missing, and extra particles. 

Recall that in the perfect honeycomb lattice, the distance to nearest neighbors is unity, while 
the distance to second nearest neighbors is \/3- In order to calculate the Defect Measure, we 
consider only a small circular region around each particle of radius (1 + \/3)/2 ~ 1.366 so that 
only three nearest neighbors, each at unit distance, should be included. Then, using the locations 
of all particles actually located within this circular region for the given lattice, the contributions 
from each of the defects is calculated, and then summed with the specified weighting, for each 
particle in the given lattice to provide a measure of the quality of the lattice as a whole. The 
weights attributed to each type of defect may be chosen according to the desired properties for the 
self-assembled lattice. These weights may emphasize correct local densities, correct alignment of 
particles, correct inter-particle distances, or any other customized weighting. 
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Displaced Particles 



•••• 



o 



Lone Particles 



ioundary Particles 



FIG. 4: To compute the Defect Measure, the number of defects in a circular neighborhood of each particle are 
weighted and summed. The types of defects used in computing the Defect Measure are displaced, missing, and extra 
particles, as well as lone and boundary particles. 



Accordingly, the Defect Measure (DM) is given by 



N 



100 

DM := > 



p=l 



ti^displaced j ^X('^pj) ' {dpj — 1-0)" 



+ '^^missing ' ^p, missing ~l~ '^^extra ' ^p, extra ~l~ ^lone ' '?p,lone ~l~ '^boundary " ^p, boundary 

where the w's are the weights attributed to each type of defect; the rip's are the integer number 
of missing and extra particles within the small circular region surrounding particle p; the r]p's are 
indicator functions that equal unity if particle p is a lone or boundary particle and zero otherwise; 
the index p ranges over all the particles in the given lattice; the index j ranges over the particles 
contained within the small circular region surrounding particle p; and dpj is the positive distance 
between particle p and particle j. As with the Template Measure, the scaling factor of 100/A^ is 
included so that the quality values are normalized by the number of particles and fall within a 
convenient range. In the first term, x(-) is a smooth cutoff function that decreases from 1 to 
at the outer edges of the circular region, or more precisely, as its argument increases from 1.275 
to 1.366. This smooth cutoff function ensures that the Defect Measure remains continuous with 
respect to motion of the particles and as the number of particles entering the circular region of 
each particle fluctuates. 

More explanation of the Defect Measure and its applications can be found in [2^, where there 
is also a discussion of other lattice quality metrics. 
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For the purposes of this paper, we use the following weight values for the Defect Measure: 



''-'displaced 1.0 , 



missmg 



1.5 , Wextra — 0.8 . 



2.0 , ^boundary — 0.1 . 



Sample values of both the Template Measure and the Defect Measure for four lattices of varying 
quality are shown in Figure [5j 




(a)TM = 131.6 
DM = 162.8 



(b)TM = 63.6 
DM = 44.2 



(c)TM = 33.9 
DM = 17.5 




(d)TM = 19.9 
DM = 12.4 



. •• 

• 

••• 

• •••••• r ~ 

•• 



(e)TM = 10.2 
DM = 4.3 



FIG. 5: Sample values of the lattice quality metrics are shown for a range of lattices. In each figure, the value of 
the Template Measure (TM) and the Defect Measure (DM) are provided. Both measures assign a lower value to a 
lattice of higher quality. 



IV. OBJECTIVE FUNCTIONS 

A natural way in which to find solutions to the self-assembly problem is to search for a potential 
that optimizes an appropriate objective function. The objective function must be well-chosen so 
that 1) optimizing the objective correlates well with the formation of the desired lattice, and 2) 
the objective function is not prohibitively expensive to evaluate. In this section, we present the 
objective functions that will be used in the potential generation methods under comparison. 
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A. The Lindemann Parameter Objective Function 



As noted in Section [III the Lindemann parameter is a measure of how much a lattice configu- 
ration has deviated from the perfect lattice configuration after a short simulation. Consequently, 
the Lindemann parameter cannot be used as a metric for lattice quality per se since it measures 
how far a lattice has deviated when initially placed in the perfect target lattice configuration, 
whereas a lattice quality metric must assess the quality of a lattice obtained after a slow cooling 
process from an arbitrary initial condition. The Lindemann parameter can nevertheless be used 
as an objective function that must be minimized in order to find potentials that lead to formation 
of the target lattice - an approach that was first proposed and developed by 2^. Intuitively, 



minimizing the Lindemann parameter produces potentials that stabilize the honeycomb lattice to 
thermal agitation. Used as an objective function, the Lindemann parameter is advantageous since 
it is much faster to evaluate than other objective functions that require much longer molecular 
dynamics simulations. 

As implemeted in this paper, the Lindemann parameter objective function is computed for a 
given interaction potential by placing 72 particles in the honeyomb lattice formation, and then 
performing a brief simulation at a temperature very near the melting temperature of the lattice. 
The use of 72 particles allows for the construction of a lattice consisting of 30 honeycomb cells that 
form an infinite honeycomb lattice when the 72-particle configuration is used to tile the plane. The 
value of the Lindemann parameter is calculated using the initial and final lattice configurations. 
The wall clock time to compute the Lindemann parameter on a single CPU is approximately 2 
seconds. 

Optimizing the Lindemann parameter is, however, an indirect method in that the quanitity 
being optimized is not the quantity that will be used to determine the final quality of the poten- 
tial. A more direct approach is to explicitly optimize the lattice quality metrics. Although the 
quality metrics require a slow cooling simulation and are consequently more expensive to evaluate, 
optimization of the quality metrics guarantees optimization of lattice quality. Indeed, an important 
observation made in this paper is that optimizing the Lindemann parameter is only moderately 
correlated with lattice quality - potentials can be found that produce a low value of the Lindemann 
parameter, yet when tested in a slow cooling simulation produce lattices of poor quality. 



13 



B. Quality Metric Objective Functions 

We also consider direct evaluation of the Template Measure and Defect Measure quality metrics 
as objective functions. To evaluate these objective functions, we start with a regular Cartesian grid 
of 64 particles with a spacing that provides a particle density equal to the density of the honeycomb 
lattice. These particles are initialized with a temperature well above the lattice melting point 
and then this temperature is slowly reduced until the particles freeze into a lattice configuration. 
The quality metrics are computed using only the final lattice configuration. Compared with the 
computation of the Lindemann parameter, these cooling simulations are expensive to carry out; a 
single call to a quality metric objective function on a single CPU takes approximately 70 seconds. 

It should be noted that the computational expense associated with evaluation of the objective 
functions arises almost entirely from the molecular dynamics simulations. Computation of the 
actual value of the Lindemann parameter, the Template measure, or the Defect measure after the 
final configuration has been obtained, requires less than a tenth of a second. 

C. Properties of the Objective Functions 

The objective functions presented have important features that influence the effectiveness of the 
optimization schemes employed. One of the chief contributions of this paper is the use of a trend 
optimization scheme that is better suited to these objective functions over a standard simulated 
annealing optimization procedure. 

The fact that the objective functions require evaluation times on the order of seconds and 
minutes implies that any optimization method that requires many thousands of evaluations to 
search the four-dimensional parameter space will necessitate a computation time measured in 
hours and days. The time taken to run an optimization algorithm will be spent almost entirely 
evaluating the objective functions, and overhead computation required by the optimization scheme 
in choosing the next location at which to evaluate the objective, for example, is practically negligible 
in comparison. The trend optimization method we propose is particularly well-suited for this 
situation in which the objective functions are expensive to evaluate. 

Figure [6] displays repeated evaluations of the Lindemann parameter objective function as each 
of the parameters in VHc(^j ooi oi) 02i os) is varied in turn, while the remaining parameters are 
held fixed at the values provided by Rechstman et al (namely, ao=5.89, ai=17.9, a2=2.49, and 
a3= 1.823). Each circle in the plots corresponds to a single evaluation of the the Lindemann 
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parameter objective function. The red line in each plot represents a trend line that is computed 
by averaging over 100 samples taken at each of 600 regularly spaced intervals along the axis of the 
varied parameter. 




FIG. 6: Each circle in the figures above corresponds to a sample evaluation of the the Lindemann parameter objective 
function. The parameter shown on the a;-axis is varied over the indicated domain, while the remaining parameters are 
held fixed at the values provided by Rechstman et al (ao=5.89, ai = 17.9, a2=2.4Q, and a3=1.823). The vertical axis 
in each plot represents the value of the Lindemann parameter. Note that the vertical scale in each case is different. 
The red line is computed by averaging over 100 trials at each location, and reveals a smooth trend in the data. 

Our first observation is that due to randomness induced by the initial conditions and simulation 
at constant temperature, the objective function does not produce repeatable values for a fixed 
set of parameter values. Repeated evaluation of the objective functions for the same parameter 
values leads to a wide range of output values. For this reason, the objective functions are not 
strictly functions, although we continue to use this term. More correctly, due to the randomness 
introduced by the initial conditions and the simulation with a thermostat, we can think of the 
objective functions as assigning a one dimensional probability distribution to each set of values 
that define the potential. Thus, a call to the objective function returns a value drawn from the 
probability distribution specified by the parameters. We shall see that the trend optimization 
method exploits the fact that the expectation values of the probability distributions vary smoothly 
with respect to changes in the parameters of the potential. 

Examining the evaluations of the objective functions reveals that they are noisy with respect 
to variations in the parameters. The functions are not smooth and certainly no derivatives of the 
objectives are available. Nevertheless, averaging over many evaluations for a fixed set of parameter 
values reveals that the objective functions do possess a smooth and slowly varying trend. As 
indicated in Figure [6l the smoothness of the average values (shown in red) indicate that, when 
averaged, the objective functions do admit a sensible notion of a minimum. 
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In summary, the salient features of the objective functions in the self-assembly problem are the 
following: 

1. They are expensive to evaluate; 

2. They are highly variable — repeated evaluations of the objective functions for the same input 
potential yields a broad range of values; 

3. They are highly non-smooth with respect to changes in the parameters; 

4. They have a smooth trend when averaged over many evaluations. 

Objectives with these features are often encountered in optimal design problems in which evalu- 
ation of the objective function for a specific set of parameter values requires completion of an actual 
laboratory experiment, or the execution of a computationally expensive simulation using random 
initial conditions 
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29|. Optimization of these objective functions is clearly not tractable using 
standard gradient-based methods. In this paper, we implement an optimization scheme ideally 
suited for objective functions of this type that has allowed us to construct an efficient procedure 
for solution of the self-assembly problem, an approach that we call trend optimization. 

V. OPTIMIZATION METHODS 

In the previous section, various objective functions have been introduced that can be used in the 
search for lattice-forming potentials. In this section, we present two methods for finding o ptim al 
values of these objectives. First, we briefly discuss simulated annealing which was used by 23|] in 
the baseline approach. Second, we present in detail the trend optimization approach. 

A. Simulated Annealing 

Simulated annealing mimics the ability of a thermal process to search a configuration space to 
find the ground state. We begin the search by evaluating the objective at an initial point, xq, in 
the parameter space and denote the value of the objective function at this location by Eq. The 
search then proceeds by jumping to new locations in a random walk. The probability, p, of making 
the jump from Xi to a new location is determined by 



p 



1 if £"1+1 < Ei , 

(5) 
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The small probability of moving to a location where the objective is in fact higher than the 
current value, allows for the search to escape from a local minimum. In analogy with true simulated 
annealing in metals, the temperature, T, is slowly reduced so that at first the search eagerly 
traverses the search space by easily escaping from local minima, and then settles to a specific 
local minimum as the temperature decreases. This simulated annealing optimization scheme was 
implemented using the GNU Scientific Library GSL_SIMAN_SDLVE( ) routine 

When applied to the self-assembly problem, the simple simulated annealing method described 
here fails to converge to a sensible minimum because of noise in the objective functions. Under- 
standably, the method gets stuck in ephemeral local minima that appear and disappear due to 
noise in the objectives. Convergence can be obtained if the objectives are sufficiently smoothed. 
Smoothing necessitates averaging over at least 20 independent simulations and thus requires added 
computational expense, and even then a large proportion of optimizations fail to converge within 
a reasonable number of objective evaluations. 

Admittedly, the simulated annealing method described here represents a very simple approach, 
and more complex methods involving adaptive search, for example, could be employed. The 
straightforward simulated annealing approach serves, therefore, as a modest but easily understood 
baseline method. We have investigated the simulated annealing method using many different 
cooling regimens including adaptive methods, and have found little improvement in each case. 

B. Trend Optimization 

Trend optimization is well-suited for problems in which the objective functions are expensive 
to evaluate, derivatives of the objectives are not available, and the objectives are noisy yet exhibit 
a simple underlying trend when averaged over many evaluations. 

In the trend approach, we use a few well-distributed evaluations of the objective function to 
generate a smooth global approximation to the averaged objective function. This smooth approx- 
imating surface, referred to as the surrogate, attempts to capture the underlying smooth trend in 
the noisy objective evaluations. The optimization then proceeds by finding the optimal values of 
this surrogate function that is computationally cheap to evaluate. Trends in the surrogate quickly 
reveal regions of the parameter space in which the optimal parameters are most likely to be found 
and hence the search is greatly accelerated. The insight provided by the smooth surrogate function 
then informs the choice of parameters at which subsequent evaluations of the objective should be 
made. 
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Booker et al |3ll ] combined the speed and facility with which the trend approach zooms in on 
regions of optimal parameter values with the rigorous convergence guarantees of patterned search 
methods developed previously by Torczon [34], to produce what is known in the literature as the 



Surrogate Management Framework. In this two-pronged approach, the trend method is used in the 
global search step of the patterned search method to accelerate the search, while the use of a polling 
step on a patterned conceptual grid provides the guarantees of convergence. In this seminal paper, 
Booker et al also applied the Surrogate Management Framework approach to the optimal design 
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of a helicopter rotor blade with thirty-one design variables. More recently, Audet et al 
have provided a generalization of Torczon's pattern search method, which they refer to as Mesh 
Adaptive Direct Search, and have applied it to optimization of the chemical treatment of discarded 
potliners to minimize release of toxic waste in the production of aluminum [35[. Mesh Adaptive 



Direct Search has since been incorporated into the Surrogate Management Framework by Marsden 
et al in [36 1. 



The range of design problems to which the trend optimization approach has been applied is 
starting to grow. Marsden et al have used the Surrogate Management Framework in the optimal 
design of airfoils to reduce noise generated in the trailing turbulent flow [s^], and a computational 



framework has been provided in 



36( 1 for optimizing design of surgeries for improved blood flow 



and cardiovascular geometry. Siah et al, have used Kriging surrogate models to design optimal 
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configurations and shapes of automobile antennae to minimize electromagnetic coupling 
Raza et al have compared methods for generating surrogate functions in a design problem seeking 
the optimal arrangement of fuel rods in a liquid metal reactor [39 1. 



A unifying theme in all these applications is the parsimonious way in which trend optimization 
is able to optimize expensive, noisy objective functions. Moreover, trend optimization is robust 
to noise in that the trends approximate the general shape of the objective function with smooth 
surfaces that quieten the noise and anomalous evaluations of the objective function. Hence, trend- 
based approaches survey the parameter landscape for large depressions, and are not distracted 
by superficial deep spikes that may arise due to noise. Because the surrogate prioritizes regions 
that consistently perform well, rather than a single "flash-in-the-pan" evaluation, the parameter 
values returned by the trend are robust to uncertainties and are more likely to reliably reproduce 
near-optimal values of the objective upon repeated evaluation. 

Trend optimization can be performed in a coarse-to-fine hierarchical manner by recursively 
building a hierarchy of trend-fitting surfaces. Each successive iteration of the procedure yields a 
new trend that focuses on the most optimal region of the search space. Successive trend surfaces 
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utilize all objective function evaluations obtained in previous iterations to more accurately model 
the objective function. After an initial global trend is developed, the search area is refined to the 
area surrounding the global minimum of the surrogate - recall that since the surrogate is smooth 
and cheap to evaluate, the global minimum of the surrogate can easily be found. Refinement of the 
search area helps to ensure that subsequent function evaluations are chosen in locations that are 
most relevant and promising. As the search area becomes more refined, successive iterations may 
use a larger basis of fitting functions, or use more sophisticated trend construction methods to more 
accurately pinpoint the location of the minimum. These features of iterative trend optimization 
yield a hierarchy of coarse to fine trends that enable the method to initially make large strides 
toward the optimal value, and then to focus ever more tightly on the exact location of the optimal 
value. 

In our discussion thus far, it remains to be described how the surrogate functions are gener- 
ated from a small number of function evaluations. This topic lies within the province of data 
approximation and fills a large body of literature. Needless to say, there are a great number of 
interpolation and fitting approaches available. Popular methods include polynomial interpolation, 
splines, Kriging, distance-based interpolation, linear and nonlinear regressions, radial basis func- 
tions, neural networks, and kernel-based approaches 



the monograph by Hastie 
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42| . For a thorough survey, please see 



Reviewing the Lindemann parameter objective evaluations depicted in Figure [6l we plainly see 
that an interpolating scheme is not appropriate for the noisy objective functions of the self-assembly 
problem. For this reason, we have chosen the ridge regression method for generating surrogate 
functions that is particularly well-suited for noisy data in high dimensions. Ridge regression was 
originally developed by Tikhonov (hence the method is sometimes referred to as Tikhonov regular- 



ization) for ill-conditioned linear regression problems [44j]. His approach was to introduce diagonal 
stabilization to regularize the linear interpolation system. The added regularization improves the 
conditioning, and introduces smoothing. 

In the current context, we are provided with a vector of M noisy measurements, y = 
[yir " tUm], of the objective function at the vector of locations x = [xi,--- ,xm]- We want 
to find a smooth function that best represents the smooth trend in this data. This trend, denoted 
T{x), is constructed as the weighted sum of basis functions: 

M 

T{x) = ^Cfc$(x,Xfc) , 

k=l 

where we must now solve for the vector of coefficients c = [ci,--- ,cm]- Gaussian radial basis 
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functions are chosen for the basis since they are well-suited for regression problems in high dimen- 
sions 



451, 
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47l |. To be specific, we choose radial basis functions of the form 

$(x,Xj) = (j){\\x - Xj\\2) 

where 

Proceeding in the standard manner for linear regression, the vector of coefficients, c, in the expres- 
sion for the trend are obtained by solving the linear system 

Ac = y (6) 

where the elements of the square symmetric matrix A are given by 

Ajj := <l>(xj,Xj) . 

This procedure assumes that A is full rank. In ridge regression, the conditioning of A is improved 
by adding diagonal regularization. The linear system in ^ above is replaced by 

A + ±l)c = y (7) 

in which u; G M, and I is the identity matrix. Regularization is obtained at the expense of introduc- 
ing the new free parameter uj. In the traditional usage of ridge regression, the analyst must choose 
an optimal value of uj that balances the need for improved conditioning with the desire to keep 
the departure from the least-squares solution small. In our context, since we are not immediately 
concerned with conditioning, we use to to control the amount of smoothing introduced. Smaller 
values of uj yield larger smoothing, while larger values of uj ensure increased pointwise accuracy to 
the noisy data. In practice, this approach is straightforward and robustly produces smooth trends 
to noisy objective functions. In Figure [71 an illustration of smooth trends generated using ridge 
regression are given for noisy evaluations of the Lindemann parameter. It must be remembered 
though that these surrogates are constructed as one-dimensional curves for illustration, whereas in 
the full self-assembly problem the surrogate functions are smooth four-dimensional hyper surf aces. 

The ridge regression method can be extended to include adaptive control of the amount of 
smoothness introduced in response to local conditions. This local ridge regression approach is 
implemented by replacing Equation ([7]) with 



A + diag 



1 1 

2uJi ' ' 2ujm 
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FIG. 7: Generating surrogate functions for noisy data. Each circle represents an evaluation of the Lindemann 
parameter as in Figure |6l The solid trend line is generated using Gaussian radial basis functions and ridge regression. 

where the vector of M free parameters [cji,--- ,00 m] can be chosen independently to adjust the 
amount of smoothing local to each measurement. In this way, a low value of uj can be chosen to 
locally increase the smoothing in a region in which the noise in the data is high, and similarly, a 
large value of oo can be chosen for measurements in regions in which there is high confidence in 
the data and noise is low. This facility is not currently implemented in our trend optimization 
algorithm, but is mentioned here to indicate the flexibility that the ridge regression method affords. 

As previously mentioned, there are many candidate approaches for trend fitting besides ridge 
regression that could be uses in the self-assembly problem. For instance, simple quadratic fitting is 
very easily implemented and guarantees convexity of the surrogate. In practice, this method also 
works remarkably well for the objectives of the self-assembly problem. 

Trend Optimization Algorithm used in the Self- Assembly Problem 

In using the trend optimization approach for the design of the potentials in the self-assembly 
problem, we have implemented it without coupling to a Direct Search method. Consequently, 
we lose rigorous guarantees of convergence, but in practice trend optimization alone performs 
remarkably well and reliably converges to an optimal value, as will be demonstrated over repeated 
trials. 

Throughout the computations, we scale the four-dimensional parameter search space to form 
a unit hypercube. We use the hierarchical approach described above with three levels of recur- 
sion. Hence, during each optimization three trend surfaces will be constructed at finer and finer 
resolution. The basic steps of the procedure are as follows. 

Let Ui denote the four-dimensional unit hypercube of the parameter space to be searched, and 
let H be the number of levels in the trend hierarchy (we use H = 'i). 
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FOR each iteration in the hierarchical trend optimization method, indexed by A: = 1, • • • ,H: 

Step 1 . Generate Sample Locations: 

Select M points, [x^f, • • • , x^j] , from Uk using Latin Hypercube Sampling 

Step 2. Evaluate the Objective: 

Evaluate the objective at each of the M locations x^, and store each corresponding 
result in y^. 

Step 3. Build the Trend: 

Use all data ^{xi,yi) : j = 1, - ■ ■ , k ; i = 1, ■ ■ ■ , m| obtained during the optimization 
so far to construct the trend surface via ridge regression. 

Step 4. Optimize the Trend: 

Quickly find xj, the location in the parameter space that globally minimizes the trend 
surface T^. 

Step 5 . Refine the Search Domain: 

Generate a new search domain, Uk+i, by reducing the size of the current search domain, 
Uk, by a factor of 2 along each dimension centered about the point x^. 

END 

After H iterations, declare x* := x^ as the parameter location that minimizes the objective 
function. If desired, the objective function can be evaluated repeatedly at x*, and then averaged, 
to generate y*, the expected value of the objective at x*. 

The most computationally expensive step is Step 2, the evaluation of the objective function. 
In comparison, optimization of the surrogate performed in Step 4 is extremely fast. To effect 
Step 4, we simply evaluated the surrogate at 5000 points randomly distributed throughout the 
search space and selected the point that produces the lowest value of the surrogate. The surrogate 
provides repeatable values so there is no need to average over many evaluations of the surrogate. 

At the completion of the algorithm, the total number of required objective function evaluations 
is M ■ H, that is, M objective evaluations during each of the H trend iterations. An important 
point is that, at each iteration, the M objective evaluations required by Step 2 are independent, 
meaning, therefore, that these M evaluations can be computed in parallel. In our optimization 
source code, we have parallelized Step 2 so that the total wall clock time required for the entire 
optimization algorithm to complete is H ■ Tobj, where Tobj is the CPU time required to run the 
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optimization with a single objective function evaluation. An important point, therefore, is that if 
the number of available computing processors is greater than M (and we consider any overhead 
communication costs between the M processors as negligible), then the execution time of the 
algorithm is independent of M. Thus, trend optimization provides an effective method to harness 
parallel computation resources for fast global optimization. 

Since in our specific implementation of the trend algorithm we recursively generate three sur- 
rogate functions, and since the evaluation of even the most expensive objective function is approx- 
imately a minute; the trend algorithm terminates after just three minutes of computation on 40 
processors having used information gathered from 120 objective evaluations. The CITerra cluster 
at Caltech has 4096 processes that could conceivably be used to find an optimal solution; however, 
we have observed that trend optimization provides reliable convergence when M (the number of 
objective evaluations made at each step) is ~ 40 or above, so that a far more modest number of 
processors is actually required. In summary, the trend optimization scheme as we have applied 
it to the self-assembly problem is able to provide parameters that optimize the most expensive 
objective functions in just over three minutes when run on a cluster of 40 parallelized processors. 

It should be noted that the hundredfold speed-up obtained by trend optimization over simulated 
annealing in the time taken to generate potentials (as stated in the introduction), is assessed 
using the total CPU time summed over all processors, and not the wall clock time. The speed- 
up is attributed to the robust manner in which trend optimization accelerates search of a noisy 
objective with a smooth trend. Hence, the speed-up attributed to the facility with which the trend 
optimization admits parallelization of the computation represents an additional time savings over 
and above the hundredfold speed-up. 

VI. METHODS FOR GENERATING POTENTIALS 

The five methods for generating potentials that we compare in this paper are described below. 
The first is a heuristic geometric method that requires no computation. The other four methods 
utilize an optimization procedure. 

A. Geometric Method 

The geometric method (abbreviated as GM) that we present here is an optimization- free pro- 
cedure that exploits differences between the geometry of the desired target lattice and competitor 
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lattices that we hope to discourage. The design of the potential is based on four main principles: 

GMl. The potential must have local minima located at each of the radial distances to nearest 
neighbors in the desired lattice and nowhere else; 

GM2. The potential may include local maxima at radial distances at which competitor lattices 
have nearest neighbors but the target lattice does not; 

GM3. When the target and competitor lattices have nearest neighbors at the same radial dis- 
tances, the energy levels of the local minima identified in GMl are chosen to energetically 
prefer the target lattice; 

GM4. A Lennard-Jones potential is spliced into the potential at the origin to provide a hard 
core potential. 

GMl and GM2 use information about the geometric structure of the target and competitor 
lattices to stabilize the target lattice, and to discourage a competitor lattice that has different 
distances to nearest neighbors. GM3 uses energetics to discriminate between lattices that have 
identical distances to nearest neighbors and differ only in the numbers of particles at those distances. 

When the target lattice is the square lattice, and the identified competitor is a triangular 
lattice, the fact that these lattices have different distances to nearest neighbors makes this method 
of potential generation a very robust approach. By explicit construction, the locations of the local 
minima and maxima discourage the triangular lattice and stabilize the square lattice. The shape of 



the potential generated for this case is shown in Figure 8(a) Notice in particular that the potential 
has a very simple form with local minima at distances of 1.0 and \/2 to encourage the square lattice, 
and a local maximum located at 1/3 to discourage the triangular lattice. In simulation this potential 
performs remarkably well. The potential robustly forms the square lattice, and defects (except for 
voids) are seldom observed. A sample square lattice obtained with this potential is shown in 
Figure |8(b)[ 



The honeycomb lattice represents a more difficult self-assembly problem since both this target 
lattice and the competitor triangular lattice have identical distances to nearest neighbors. The 
only difference between these lattices that must be exploited is the different number of neighbors 
at each distance. Hence, in GM3, the relative heights of the local minima are chosen to ensure 
that the target lattice represents a lower energy state than the competitor lattice and is thus more 
likely to form upon cooling. 
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FIG. 8: Potentials developed using the heuristic geometric method for both the square and honeycomb lattices are 
shown in Figure (a). Figures (b) and (c) show typical lattices obtained using these potentials. 

The application of the GM method for the specific case of constructing a potential that favors 
the honeycomb lattice and discourages the triangular lattice proceeds as follows: 

1. Recognize that both the triangular and honeycomb lattices have nearest neighbors at a 
distance of 1.0, and second nearest neighbors at a distance of \/3. The number of neighbors 
at those distances (the lattice coordination numbers) for the honeycomb lattice are 3 and 6, 
while the triangular lattice has 6 and 6. 

2. Assign local minima at distances of 1.0 and Vs. 

3. Since the triangular lattice has more particles at distance 1.0 than the target lattice, raise 
the first minimum at 1.0 with respect to the height of the minimum at -v/3- 

4. Splice in a 6-12 Lennard-Jones type repulsive potential at the origin to simulate a rigid core. 

5. Use cubic splines to piece together the repulsive core at the origin and the local minima at 
their various locations and heights. 

One other consideration that we include when designing this potential is that we are careful 
not to make the local maximum between the first and second minima too high. This ensures that 
particles that are stuck in the triangular formation have increased likelihood of escaping to the 
lower potential well of the honeycomb lattice. 

The cubic splines are constructed so as to have zero derivative at the local minima, thus ensuring 
that the potential is continuously differentiable on the whole positive real line. A sample potential 
constructed using this method for the generation of honeycomb lattice potentials is shown in 
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Figure 8(a) , Notice that the potential has local minima located at distances of 1.0 and a/S, with the 
first minimum located above the second minimum. Since the triangular lattice has more particles 
located at distance 1.0 than the honeycomb lattice, this ensures that the triangular lattice is a 
higher energy state than the honeycomb lattice. A sample lattice produced using this potential is 



shown in Figure 8(c) , The fact that this heuristic geometric method is able to produce honeycomb 
lattices of this quality without having to use any computation and optimization is quite remarkable 
and until now has been overlooked in the self-assembly literature. 

We refer to the geometric method as a heuristic method since we do not provide formulas or 
algorithms for choosing the relative heights of the local minima and maxima. The only prescribed 
constraints are the locations of the local minima that, by construction, ensure stability of the target 
lattice. In practice, we find that any sensible choice of the heights that avoids sharp gradients in 
the potential works reasonably well. In the implementation of our computer code, the user need 
only input the location and heights of the desired minima and maxima, and the spline fitting and 
splicing of the Lennard-Jones potential are computed automatically, which makes the approach 
very simple to execute. Given the simplicity of the geometric method, the ease with which it 
is implemented, and the relatively high quality of the resulting lattices, it represents a "rough- 
and-ready" approach that anyone needing to design a potential would do well to consider before 
pursuing more computationally expensive schemes. 

Whereas the geometric method utilizes only the geometries of the static target lattice and 
competitor lattice, the optimization methods discussed next incorporate information about particle 
dynamics by optimizing parameters in the potential with respect to dynamic particle simulations. 

B. Baseline Method 



23( 1 uses the simulated annealing optimization 



The baseline simulated annealing method of 
procedure with the Lindemann parameter as the objective function. As mentioned previously, the 
simulated annealing procedure fails to converge to a meaningful minimum unless the objective 
function is averaged over many trials. We refer to this method with the abbreviation SA-LP20 to 
indicate that the simulated annealing approach is applied to the Lindemann parameter averaged 
over 20 independent evaluations. The large number of samples required to sufficiently reduce 
the noise in the objective function, means that it is impractical to apply the simulated annealing 
optimization procedure to the more expensive quality metric objective functions. The simulated 
annealing procedure is initiated with initial guess solutions chosen randomly and uniformly from 
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Summary of Potential Generation Methods 


Method 


Optimization 


Objective 


Cost (s) 


GM 


Geometric method 





SA-LP20 


Sim. anneal. 


Lindcmann (20 sample avg.) 


40 


T-LP 


Trend 


Lindcmann 


2 


T-TM 


Trend 


Template Measure 


70 


T-DM 


Trend 


Defect Measure 


70 



TABLE I: Each of the five methods for generating potentials is listed here with the associated computational cost 
of evaluating the objective function. Notice that the Geometric Method is not an optimization-based method and 
incurs no computational cost. 

the parameter space described by the inequahties in hne ^ of Section [Til 
C. Trend Optimization Methods 

The remaining three methods that we consider utihze the trend optimization method apphed to 
the Lindcmann parameter, the Template Measure, and the Defect Measure as objective functions. 
These methods are abbreviated T-LP, T-TM, and T-DM respectively. Since the trend optimization 
method is well-suited to noisy objective functions, the objectives do not need to be averaged over 
many runs when evaluated. 

A summary of all the potential generation methods under study is provided in Table HI 
VII. RESULTS 

In this section, we compare the effectiveness of the proposed methods for generating interaction 
potentials that lead to the self-assembly of the honeycomb lattice. In particular, in accordance 
with the comparison criteria enumerated in Section [Til we seek answers to the following questions: 

1. How much computational effort is required to generate the potentials? 

2. What is the quality of the lattices generated by the potentials? 

3. How reliably do the potentials form high quality lattices given uncertainty in the initial 
conditions of the particles? 

Answers to these questions, as well as a comparison of the five potential generation methods 
discussed in Section IVH are summarized in the suite of plots shown in Figures M and [lOl In 
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brief, Figure [9] shows the superior effectiveness of trend optimization over simulated annealing in 
minimizing the respective objectives as a function of the number of objective evaluations required, 
while Figure shows the quality of the lattices obtained by the potentials that were generated by 
the optimizations. 




100 1000 10000 

Number of Function Evaluations 



(a) Simulated annealing optimization of 
the Lindemann parameter (SA-LP20). 



50 

45 

40 

S 35 

I 30 
a> 

o 25 

Q. 

^ 20 
15 
10 
5 
0, 




Number of Function Evaluations 



(c)Trend optimization of the Template 
Measure (T-TM). 
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(b)Trend optimization of the 
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(d)Trend optimization of the Defect 
Measure (T-DM). 



FIG. 9: Averaged objective function values versus number of objective function evaluations for four optimization 
methods are shown here. In the results for simulated annealing shown in (a), dashed lines indicate optimizations 
that failed to find an optimal value. The single bold line represents the mean of the optimizations that did converge 
to an optimal value (indicated by solid lines). Most notably, trend optimization reliably converges to a minimum 
value of the Lindemann parameter with a one- hundredfold reduction in the number of required objective function 
evaluations when compared with the simulated annealing approach. 



The geometric method (GM) is the simplest method to implement. It does not require any 
objective optimization to search for parameters and consequently does not incur any computational 
cost. In order to evaluate the quality of the lattices produced by the geometric method, one-hundred 



cooling simulations were performed using the honeycomb potential shown in Figure 8(a) After 



each simulation, the Template Measure and the Defect Measure were computed. Averaging over 
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all these runs yielded the following scores: 

1. Template Measure = 17.9 with standard deviation 3.4, 

2. Defect measure = 11.8 with standard deviation 4.5. 



Comparison with the lattice qualities obtained using the optimization-based methods as shown in 
Figure [TOl indicates that these quality scores are exceptionally good, and that assiduous computa- 
tional optimization is required to produce potentials that produce lattices of higher quality. 

The four optimization-based methods (SA-LP20, T-LP, T-TM, and T-DM) were each run for an 
increasing number of allowed objective function evaluations. For each number of objective func- 
tion evaluations, each trend optimization method was executed in twenty independent trials, with 
each trial generating parameters for the interaction potential that seek to minimize the respective 
objective function. In Figure [H the average of the values of the objective function obtained in 
the twenty trials is plotted for each number of function evaluations (the error bars indicate the 
standard deviation over the twenty trials for each number of function evaluations). 

After all the optimizations are completed and the potentials have been generated, the quality 
of each potential must be tested. For each number of function evaluations, and for each of the 
twenty independent trials, the potential produced by a method was quality tested by running 
twenty cooling simulations on a system of 225 particles, and then measuring the quality of the final 
lattices using both the Template Measure and the Defect Measure. For a single trend optimization 
method, this requires 

f 12 different numbers of function evaluations ^ / 20 independent trials \ ^20 cooling simulations ^ 



each trend method / V each number of function evaluations / V each independent trial / 

4800 cooling simulations \ 



each trend method / 



In other words, 240 independent optimization trials are required to produce a single curve in Figure 
[9l and 4800 cooling simulations are required to produce a single curve in Figure [101 The cooling 
simulations were initialized with a temperature approximately 1.5 times the melting temperature 
of the lattice and then slowly cooled using a Nose-Hoover thermostat to less than ten percent of the 
melting temperature. At the completion of each cooling simulation, both the Template Measure 
and the Defect Measure were computed to measure the quality of the resulting lattice. Averaging 
over the 20 cooling simulations yields two quality scores for each potential - one for each quality 



metric. Results for the Template Measure are shown in Figure 10(c) while results for the Defect 



Measure are shown in Figure 10(d) 
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After reviewing the results in Figure O we see that the trend approach clearly provides a faster 



and more robust method for optimizing the objective functions. Figure 9(a) shows the results 
from twenty simulated annealing optimizations for different randomly chosen starting points in 
the search domain, and indicates that the success of the method is highly variable. Simulated 
annealing often fails to find optimal parameter values and remains stuck in local minima even 
after 10,000 evaluations of the objective. Of the twenty simulated annealing optimization trials 
performed, only nine were able to find an optimal value of the Lindemann parameter less than 



0.15. The bold line in Figure 9(a) represents the mean over only these nine best-performing 
trials. Moreover, many of these simulated annealing optimizations that do converge require more 
than 6,000 objective evaluations. In contrast, we see from Figure |9(b)| that trend optimization 
reliably finds optimal values after sixty evaluations of the Lindemann parameter. We conclude 
that when trend optimization is used to optimize the Lindemann parameter, we obtain a one- 
hundred-fold reduction in computation time over the simulated annealing method of [23|], and that 



the optimization is more robust. The speed-up can be attributed to the fact that the objective 
function is noisy yet has a simple trend - two properties for which trend optimization is ideally 
suited. 

Furthermore, the accelerated search of the trend method makes it possible to use trend opti- 
mization with the more expensive quality metrics as objective functions. Doing so provides the 
extra guarantee that the generated potentials produce high-quality lattices. As indicated in Fig- 
9(c) and |9(d)"| trend optimization reliably finds optimal values of the lattice quality objectives 
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in less than 120 function evaluations, although it must be remembered that these objectives are 
approximately 35 times more expensive to evaluate than the Lindemann parameter. Nevertheless, 
the total CPU time required to perform these optimizations is still less than the time taken by 
simulated annealing to optimize the Lindemann parameter. 

Figure [To] shows the quality of the lattices generated by the potentials produced by the various 
methods, as measured using both quality metrics. Several important observations are to be made 
after reviewing these plots. First, optimization of the Lindemann parameter leads to lattices of 
modest and unreliable quality, indicating that the Lindemann parameter and lattice quality are 
only moderately correlated. Second, even when simulated annealing is successful in optimizing the 
Lindemann parameter, the corresponding potential may yield lattices of poor quality. This occurs 
because simulated annealing may find minima corresponding to very narrow wells that do not pro- 
vide robustness against uncertainty in initial conditions. In contrast, since trend optimization seeks 
out the general trend over the entire search space, this method finds minima that more robustly 
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(a) Quality of lattices generated using 
simulated annealing measured with the 
Template Measure. 
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(c) Comparison of lattice quality using the 
Template Measure. 
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(b) Quality of lattices generated using 
simulated annealing measured with the 
Defect Measure. 
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(d) Comparison of lattice quality using the 
Defect Measure. 



FIG. 10: Here we compare the quality of the lattices produced using the four optimization methods. In (a) and 
(b), the quality of the lattices produced using the simulated optimization method is shown. It should be noted that 



only the nine potentials generated by the convergent simulated annealing trials in Figure 9(a) were used to produce 
these plots. Figures (c) and (d) show the quality of lattices produced using the three trend optimization methods. 
Using trend optimization directly on the quality metrics reliably produces high quality lattices. Optimization of the 
Lindemann parameter is only moderately correlated with improved lattice quality. In these Figures, recall that values 
of the Template Measure and the Defect Measure for purely random configurations are 132 and 161, respectively, 
and that the computation-free Geometric Method produces lattice quality scores of 18 and 12, respectively. 



lead to high quahty lattices. Finally, the two lattice quality metrics are reasonably consistent in 
that a potential generated by optimizing one of the quality measures is also considered high-quality 
in the other measure. 

Figure [11] shows the range of potentials generated by the trend optimization methods. The 
potential provided by [2^ is also shown for reference. Most striking is that trend optimization un- 
covers an entirely new family of potentials not previously considered. In this family of potentials, 
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FIG. 11: Potentials for the generation of honeycomb lattices. Vr is the potential previously provided by Q. Trend 
optimization generates potentials similar to [231 (labelled Va), but also uncovers a new family of solutions that have 
a more repulsive shape and no local minima (sample potentials in this family are labelled Vb and Vc)- Remarkably, 
these repulsive potentials robustly form large regions of honeycomb lattice without defects. 



the exponential term is dominant and consequently, their shape does not admit a local minimum. 
The repulsive shape of these potentials leads to higher quality lattices since particles do not get 
stuck in local minima associated with local potential wells. In simulation, we observe that the de- 
fects continue to move through the configuration until they leave the domain entirely, or annihilate 
one another through collisions. 

Figure [12] shows two sample lattice configurations obtained for 4128 particles simulated using 
the potential labelled Vc in Figure [TT] that was generated with trend optimization on the Tem- 
plate Measure. The parameter values for this potential are ao=5.771, ai=23.594, 02=0.574, and 
03=1.816. The final lattices exhibit large regions of extremely well-formed honeycomb lattice that 
we have not previously observed using potentials that contain local minima. In Figure 12(a), a 



prominent grain boundary lies between two large areas of well-formed honeycomb lattice. As in 

iso- 



Nature, this grain boundary forms when the cooling is not sufficiently slow. In Figure 12(b) 



lated defects are visible in the lattice; however, it should be noted that these defects are not yet 
"frozen" into the lattice, and given enough time will eventually leave the domain or disappear 
through mutual collision. 

The repulsive potentials found using trend optimization are remarkably effective at producing 
large regions of almost defect free honeycomb lattice despite their relatively simple shape. The 
effectiveness of these potentials suggests that a far simpler basis of potential functions can be 
used in the parameterization of Vhc(^)- In these potentials it is the exponential decay term that 
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(a) (b) 

FIG. 12: Honeycomb lattices formed using the repulsive family of potentials discovered using trend optimization, 
(a) Two large regions of extremely well-formed honeycomb lattice meet to form a grain boundary, (b) Large regions 
of the honeycomb lattice are formed with a few isolated defects. These defects continue to move, even after the 
lattice is frozen, until they exit the boundary or are annihilated by collisions with other defects. 



dominates over the Gaussian term. The numerical simulations seem to suggest, and it would be 
interesting to pursue rigorously, that for a fixed density, the honeycomb lattice is a global minimizer 
(or ground state) of the exponential decay pairwise interaction potential. This path is consistent 
with the numerical findings of Jagla et al in [5^ for a ramp type potential. 



VIII. ANISOTROPIC POTENTIALS 

The methods presented thus far can be used to quickly generate isotropic potentials that produce 
high quality lattices. However, the potentials for self-assembling honeycomb lattices are not robust 
to variations in density - the lattices only form if the initial density of the particles is very near the 
ideal density of the target lattice. If the initial density is much less than that of the target lattice, 
then a triangular lattice will be favored. 



23l ] addressed this issue by searching for optimal potentials 
over a range of densities. However, the tolerance for variation in density is still very narrow and 
may indeed be a fundamental limitation of isotropic potentials. Allowing the potentials to have 
angular dependence leads to the robust formation of high-quality lattices from initial conditions 
with large variations in density. Moreover, admitting potentials with angular dependence allows for 
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the construction of potentials that form more exotic lattices, such as the kagome lattice, which has 
not yet been accomplished with purely isotropic potentials [23(]. Certainly, the use of anisotropic 



potentials is well-motivated by abundant natural examples of anisotropic interaction potentials in 
Nature - the water molecule serving as an ubiquitous prototype. 

In order to introduce potentials with angular dependence, we must first extend the config- 
uration space of the particle system to include an angular coordinate. We no longer consider 
the particles as point particles, but rather as two-dimensional sliding disks, each with radius R 
and uniform mass density. The configuration manifold of each particle is now x S^. The 
configuration of the i^^ particle is described by the coordinate chart {xi,yi,9i) as indicated in Fig- 
ure [T3l and we write the vector of coordinates describing the configuration of all N particles as 
X = [(xi, yi, ^i), • • • , {xN,yN, On)] ■ The total interaction potential over all particles now has the 




FIG. 13: To implement an anisotropic potential, particles are now modeled as two-dimensional disks whose con- 
figuration is described by a location for the center of the disk, {x,y), as well as a heading angle 8. The interaction 
potential between two particles is a function of the distance between the particles, as well Eis the angle, a, that lies 
between the heading angle and the bearing toward the second particle. 



form 

N 

l/(x) = ^ Vpi,ir{xi - Xj,yi - yj,ei, 9j) (8) 

i<j 

where Ipair : M x M x S""*^ x 5^ — > M is the pairwise interaction potential between particles that 
depends not only on the relative displacement between particles, but also the angular displacement 
of each particle relative to the angular bearing of the other particle. By construction, we make 
Vpair symmetric with respect to particle interchange. Specifically, 

^pair(a^i - Xj,yi - yj,6i, 9j) := ip{xi - Xj,yi - yj,9i) + ipixj - Xi, yj - yi, Oj) . (9) 
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It remains to choose the functional dependence of to reflect the desired symmetry in the 

target lattice. Consider the interaction potential 

mx,Ay,9,) := ^ - ^ [l - z^sin^ (^)] (10) 

where 



r := v/(Ax)2 + (Ay)2 (11) 

is the radial distance between the particles, and 

o^k ■= ^fc ~ arctan(Ay, Ax) (12) 

is the angle between the heading of particle k and the bearing toward the second particle (see Figure 
[T3l ) This potential is simply a Lennard- Jones potential with added amplitude modulation in the 
azimuthal direction on the attractive term. The periodicity of the trigonometric functions ensures 
that the potential has n-fold radial symmetry, and that each particle has preferred directions along 
which it feels the attractive pull of neighboring particles. The free parameter z/ is chosen to adjust 
the shape of the potential. When is zero, the potential collapses to the isotropic Lennard-Jones 
potential. For values of v between zero and unity, the potential has n potential wells symmetrically 
distributed in the azimuthal direction. Taking values of v greater than unity raises the repulsive 
regions between the potential wells. Hence, v is a parameter that determines how strongly the 
anisotropic potential prefers the binding site directions. In simulations, we have observed that 
taking larger values of v produces lattices with less defects since particles that do not align with 
the preferred binding directions fall into these more repulsive regions between the wells and create 
a configuration with much higher energy. These defects are quickly removed by vibrations in the 
lattice. In practice, a compromise must be met since the larger regions of repulsion created by 
higher values of v increase the time taken for self-assembly to occur - particles only bind with 
one another if they approach each other along an ever more narrower binding direction. In the 
simulations that follow, we have used a value of = 1.5. 

By changing the integer value of n, we can induce the formation of lattices with desired n-fold 
symmetry. Surface plots of potentials with 3-fold and 4-fold symmetry (n=3 and n=4) as a function 



of the radial coordinate r and the azimuthal angle a, are provided in Figures 14(a) and 14(b) , A 



honeycomb lattice and a square lattice produced with these potentials using particles initialized at 



low density are shown in Figures 14(c) and 14(d) 



For the creation of more exotic lattices, we alter the azimuthal modulation of the potential 
even further. Each particle in the kagome lattice, for example, has binding sites at angles of 0°, 
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FIG. 14: The potential with three- fold symmetry shown in (a) yields the honeycomb lattice shown in (c). Similarly, 
the potential with four- fold symmetry shown in (b) yields the square lattice shown in (d) . 



60°, 180°, and 240° [49|. To ensure preferred binding along these directions, we must modify the 
azimuthal dependence of the Lennard-Jones potential accordingly. Before doing so, we first express 
the potential function -0 of line (jlOp in the following equivalent way: 

1 1 



ilj{Ax,Ay,9k) 



2^12 j,6 



[l-uS{ak)] 



where 



S[ak) := sm j . 



(13) 



(14) 



As written here, ip produces a potential with n-fold symmetry. In order to produce a potential 
that favors the kagome lattice, we must simply introduce a new definition for the functional form 
of5(-). 

Let B := [bi, • • • ,bn] denote the ordered list of n desired binding directions measured in radians 
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satisfying 



= 61 < • • • < 6„ < 27r . 



(15) 



Note that without loss of generahty, we may prescribe that the first binding direction lies along 
the ray corresponding to zero radians. For the honeycomb lattice, i?HC •= [O; 3"] 1 while for the 
kagome lattice we have i?kagome '■= [O, tt, ^] . Then, for a G [0, 27r), we use the elements in the 
list B to define S{a) piecewise as follows: 



S{a) :-- 



sm 



sm 



ce-bj 
bi+i~bi ' 



) if 6j < a < , for i = 1, 



2TT~b„ 



(16) 



if 6„ < a < 27r . 



When S{a) defined in this way is substituted into the expression for the anisotropic potential 
function in line ()13p . it provides the necessary azimuthal modulation of the Lennard- Jones 
potential to produce potential wells along the binding directions specified in the list S. A plot of 



S{a) using -Bkagome is provided in Figure 15(a) Notice that S{a) has local minima precisely at the 
binding site angles of 0°, 60°, 180°, and 240° consistent with the kagome lattice (recall from line 
(jlSP that minima in S{a) lead to minima in the interaction potential). Consequently, the angular 
dependence in the resulting interaction potential favors the bond structure peculiar to the kagome 



lattice. The potential and a lattice resulting from this potential are shown in Figures 15(b) and 



15(c) respectively. 



IX. CONCLUSIONS AND FUTURE WORK 



We have presented and compared methods for generating potentials that lead to the self- 
assembly of specified target lattices. In particular, we have addressed the problem of designing 
pairwise interaction potentials that induce the formation of a honeycomb lattice when a planar 
system of particles is cooled. 

We have demonstrated that reasonably high quality lattices can be produced using a heuristic 
computation- free geometric method. The geometric method provides principles that utilize purely 
geometric information to design potentials that by construction favor and stabilize the target 
lattice. 

A trend optimization algorithm has been introduced that quickly and robustly finds optimal 
shapes of the interaction potential that lead to the self-assembly of lattices of high quality. The 
success of the trend method lies in its ability to quickly locate minima in a noisy and expensive 
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FIG. 15: The function in (a) is used to azimuthally modulate the amplitude of a Lennard- Jones potential to 
encourage binding directions that favor the kagome lattice. The resulting potential surface is shown in (b), while a 
kagome lattice formed with this potential is shown in (c). 



objective function. Moreover, the potentials discovered using the trend optimization procedure 
robustly form high quality lattices with respect to variations in the initial conditions of the particles. 
We have seen that the trend optimization method has discovered a family of potentials characterized 
by very simple exponential decay profiles that routinely lead to the formation of the honeycomb 
lattice. 

Our trend optimization algorithm robustly and routinely finds optimal values of the objective, 
although it must be noted that as currently implemented, the algorithm does not provide rigorous 
guarantees of convergence. Convergence can be guaranteed, however, by incorporating a polling 
step as required in the Surrogate Management Framework. It would be interesting to investigate 
the effects of polling on the optimal results and the efficiency of the method. One expects that the 
addition of a polling step to guarantee convergence will represent only a marginal increase in the 
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overall computational cost of the algorithm. 

The geometric method has also been extended to the design of anisotropic potentials. Azimuthal 
dependence of the interaction potential allows for the formation of the kagome lattice which has not 
previously been performed using isotropic potentials. Incorporating anisotropy into the potentials 
allows for the formation of lattices over a wide range of particle densities. 

An auxiliary contribution of this paper is the development of two metrics for objective analysis 
of lattice quality. The development of these metrics was necessary for comparison of the lattices 
produced by the proposed potentials. 

We anticipate that the methods presented here will naturally extend to three dimensions without 
impediment. Rechtsman et al have recently investigated the design of potentials for self-assembly 
of three-dimensional lattice structures [51[ . Most notably, they have demonstrated the formation of 
the diamond and wurtzite lattices. As a matter of course, we intend to apply the methods presented 
here to the design of potentials in the three-dimensional self-assembly problem and expect the trend 
optimization method to provide a significant computational savings in the design of potentials. Of 
particular interest is to determine if trend optimization recovers the three-dimensional result of [511] , 



or if a new, simpler, and perhaps more robust family of solutions is discovered. 

A natural extension of the methods presented is to investigate the use of multi-specie and multi- 
body potentials for the self-assembly of quasicrystals and Penrose tilings. Generating lattices with 
prescribed geometry and structure is highly motivated by the desire to produce photonic crystals 
with specified optical properties. Also, we intend to explore the use of the trend optimization 
method for designing potentials that lead to the the formation of hexahedral meshes over compli- 
cated three-dimensional volumes. Producing hexahedral meshes is a notoriously difficult problem 
that is currently a major stumbling block for continuum mechanics computations that use a finite 
element method. 

More broadly, we expect to use trend optimization to understand more deeply the fundamental 
limitations and extraordinary possibilities of self-assembly. It will be interesting to investigate, 
for example, the role of optimal potential design in natural systems. In this regard, a study of 
self-assembly is a study of life itself. Are there mechanisms at the level of local interactions that 
allow for the differentiation of cells that self-assemble to form bone, for instance, from those that 
form brain tissue? How much control authority over the superstructure formed via self-assembly is 
provided by control over the local interactions? Are there natural systems in which a small amount 
of flexibility in the properties of the local interactions allows for large and beneficial changes in 
structure at the macroscopic level? We believe that insights to these questions may be afforded by 
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finding optimal potentials through the method of trend optimization. 
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